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In this paper we explore one of the most important features of the Galerkin method, which is to 
achieve high accuracy with a relatively modest computational effort, in the dynamics of Robinson- 
Trautman spacetimes. 
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I. INTRODUCTION 

The General Theory of Relativity (GR) consti- 
tutes one of the major scientific achievements of the 
past hundred years. Its theoretical framework es- 
tablishes that the gravitational field and the geom- 
etry of the spacetime are the same physical entity 
and that the equations ruling the gravitational dy- 
namics are field equations of the spacetime geome- 
try. GR has so far passed extremely well all exper- 
imental tests[l|, in opposition to other relativistic 
extensions as scalar-tensor theories. Basically, GR 
establishes a system of ten nonlinear, coupled par- 
tial differential equations that governs the dynamics 
of the gravitational field represented by the symmet- 
ric second order metric tensor g a p, a, /3 — 0,1,2,3. 
The nonlinear nature of the field equations is the 
main obstacle to obtain exact solutions of the field 
equations, unless we make appeal to symmetries in 
the solutions, or alternatively the use of perturba- 
tive approaches are adopted. Therefore, in order 
to study the dynamics of the gravitational field in 
more general situations the use of numerical tech- 
niques seems to be the only possible strategy to cir- 
cumvent the difficulty posed by the nonlinearities of 
the field equations. In this context numerical rel- 
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ativity has become a very fertile and at the same 
time challenging field of research as recently consid- 
ered in several interesting reviews [2j. There the im- 
provement of specific numerical techniques adapted 
to relativistic problems, along with the increase of 
computational resources figure as promising factors 
for the advance of numerical relativity. Nonetheless, 
despite all the progress made so far the complete 
understanding of important problems in relativistic 
astrophysics such as non-spherical collapse and non- 
linear regimes of emission of gravitational waves is 
still not complete. 

A promising approach to treat numerically non- 
linear problems is provided by the so-called spec- 
tral methods The spectral methods adopt a dis- 
tinct strategy if compared with the finite difference 
scheme. For instance, considering a function u(t, x) 
satisfying a given one dimensional partial differen- 
tial equation, it will be approximate as a series of 
the type u a (t,x) = J2k=o a k(t)"4>k(t), where the ba- 
sis or trial functions ipk(x) are known analytical 
polynomials such as Fourier, Legendre, Chebyshev, 
etc. In general, by increasing the truncation order 
N, u a (t,x) approaches of the exact solution of the 
problem in the mean. There are distinct types of 
spectral methods among which we list the Galerkin 
method[l,Q, the collocation method^] and the Tau- 
method. These methods have an attractive feature 
which is to transform any partial differential equa- 
tion into a finite set of coupled ordinary differential 
equations, or simply a dynamical system whose di- 
mension is dictated by the truncation order N. An- 
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other important and robust feature is the high ac- 
curacy achieved with a small truncation order, what 
corresponds to use moderate or low computational 
resources. 

In this paper our objective is to present a con- 
sistent low dimensional dynamical system approach 
provided by the Galerkin method applied to the 
problem of non-spherical collapse with the emission 
of gravitational waves, in the realm of Robinson- 
Trautman geometries [7| . In this way, the paper is 
divided as follows. The basic equations that gov- 
ern the dynamics of Robinson- Trautman spacetimes 
are presented in Section 2. In Section 3 we discuss 
and implement efficiently the Galerkin method. In 
Section 4 numerical experiments are performed in 
order to test the efficiency and convergence of the 
method. Section 5 is devoted to the conclusions and 
final remarks. 



II. THE DYNAMICS OF 
ROBINSON-TRAUTMAN SPACETIMES: 
THE NONLINEAR DIFFERENTIAL 
EQUATION 



the evolution of the gravitational field, in other 
words, it allows to evolve the initial data Kg(9) = 
K{u = uq,9) prescribed on a given null surface 
u = uq = constant (except in the case m = 0). 
Despite its relatively simple form no general exact 
non-stationary solution of RT equation is known. 
However by restricting to stationary configurations 
two important solutions are known. The first is the 
Schwarzschild solution that represents the gravita- 
tional field a black hole characterized by K = Kq = 
constant, which implies A = K^ 2 with a correspond- 
ing Schwarzschild mass M — m^K^. The second so- 
lution represents a boosted black hole[9( along the 
axis of symmetry such that 



cosh 7 + cos 9 sinh 7 ' 

where 7 is a parameter associated with the velocity 
of the boost, v = tanh7. The constraint equation 
([2|) in this case yields A = 1. 

III. THE SPECTRAL METHOD 
APPROACH 



Robinson- Trautman (RT) spacetimes are the sim- 
plest axisymmetric spacetimes that can be viewed as 
the exterior geometry of an isolated source emitting 
gravitational waves 7, 8] or the exterior of a distorted 
black hole. In a suitable coordinate system the line 
element for RT geometries is expressed as 



ds 2 = X(u, 9) 



2niQ 



r 2 K 2 (u,9)(d9 2 



„ K{u,9) \ l2 „ , , 
+" 2r ^ du 2 + 2dudr 

K(u, 9) J 

sin 2 9dip 2 ), 



(1) 



where too is a constant related to the total mass of 
the system and dot means derivative with respect to 
u. According to Einstein field equations the function 
X(u,9) is related to K = K(u,0) by the constraint 
equation 



X(u,9) = — - — + cot0, (2) 

where the subscript 9 denotes derivative with re- 
spect to the angle 9. The remaining vacuum Einstein 
equations impose that the function K(u, 9) satisfies 
the evolution equation 



6m o-| 



(Xg sin 9)g 



= 0, 



(3) 



2K 2 sin 9 

which henceforth will be denoted as the Robinson- 
Trautman (RT) equation. This equation governs 



It will be useful to express K (it, 9) as 



K(u,x) = A ei Ql - u ' x \ 



(5) 



where Aq is a constant, and for convenience we have 
introduced the variable x — cos 9 with — 1 < x < 1. 
The first step to apply the Galerkin method is to 
express Q(u,x) as a power series with respect to a 
suitable set of basis functions chosen as the Legendre 
polynomials due to the symmetry of the problem and 
the required regularity of K(u, x) with respect to x. 
In this way an approximated expression for Q(u,x) 
is established and given by 



}a{u,x) 



N 

E 

k=Q 



b k (u)P k (x), 



(6) 



where ./V is the truncation order, bk{u) and Pk(x) are 
the modal coefficients and the Legendre polynomials 
of kth order. This basis functions define an abstract 
projection space for which an internal product can 
be defined, in this 



(P 3 (x),P k (x)) 



Pj(x)P k {x)dx 



2S 



kj 



2k + 1' 



The next step is to substitute the Galerkin decom- 
position for Q(u,x) into Eq. ([2]) to obtain an ap- 
proximated expression for X(u,x), or 
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Xa(u,x) 



-Qa(u,x) 



N 



k=0 



The so called residual equation follows after intro- 
ducing the approximated expressions for Q(u, x) and 
X(u, x) given above into the RT equation ([3|). After 
some direct calculation we obtain schematically 



• The second strategy consists, on the hand, in 
applying the collocation method @ to express 
e -Q a (u,x) ag a ser i es f another basis functions, 
or 



-Q a (lt,x) 



M 

E 

fe=0 



q k {u)T k {x), 



(12) 



N 

Res(w, x) = y^bk{u)P k {x) 

fc=0 



a -Q (u,i) 

6m Al 



[(1-x 2 )A;]',(9) 



where prime stands for derivative with respect to 
x. In fact the residual equation can be a good mea- 
sure of the error of truncation and a valuable test [13] 
for the convergence of the Galerkin decomposition. 
According to the traditional Galerkin method the 
projection of the residual equation into each basis 
function P n (x), n = 0, 1, N is assumed to be zero 
in order to minimizing the error of truncation, or 
(Res(u, x), P n {x)) — 0, rendering the following dy- 
namical system for the modal coefficients 



b n {u) 



(2n- 



1) f 1 ( N 
p- _/ i cxp ( -^J>k{u)P k {x) 



12m Al 
[(l-x 2 )X' a ]'P n (x)dx 



(10) 



with n = 0, 1, ...,N. 

Here some remarks are necessary. In order to ex- 
press the rhs of the above equation in terms of the 
modal coefficients, it is necessary to integrate a large 
expression involving an exponential of a series of 
Legendre functions, the size of which depends on the 
truncation order N. This integration must be neces- 
sarily done symbolically (in Maple or Mathcmatica) , 
but unless N < 2 - a very low truncation order - 
such an integration cannot be performed. Thus, we 
can adopt two distinct strategies based on further 
approximation schemes: 

• The first is to assume 



where T k (x) are the Chebyshev polynomials, 
and qk{u) the new modal coefficients given 
as functions of the original coefficients bj(u). 
These relations are fixed by the fact that 
at the collocation points Xj = cos(irj/M), 
j = 0, 1, M the above decomposition agrees 
with the exact expression, or 



-Q a (u,Xj) 



M 

E 

k=0 



q k (u)T k (xj), 



and together with the discrete orthogonal- 
ity relation Sj k = j^J2rLo i^ T j( x n)Tk(x n ), 
where cq — cm — 2 and c k = 1, for 1 < k < 
M — 1. As a consequence each modal coeffi- 
cient q k can be expressed in terms of b k . In 
both schemes a N + 1 dimensional dynamical 
system for the modal coefficients is obtained 
after integrating the rhs of Eq. (|T0|) . 



IV. NUMERICAL EXPERIMENTS 

Our first numerical task will be to study the accu- 
racy and convergence of the approximative schemes 
outlined at the end of the last Section. To this 
end we consider the exact expressions for some ini- 
tial data Kq(x) — K(uo,x) identified as an oblate 
spheroid and a black hole in motion given, respec- 
tively, by 



E 



(-QaV 



(11) 



which renders a polynomial expression for the 
exponential. Introducing this expression into 
Eq ([9]) the integration can be done promptly 
without much cost of CPU time for reason- 
able choices of N and J. As a matter of fact, 
this additional approximation can be accept- 
able for not large values of \Q\ as we are going 
to see later. 



K (x) 



K {x) 



fcov 7 ! — 



cosh 7 + x sinh 7 ' 



(13) 



where fco and 7 are constants and e is the eccen- 
tricity. As a further assumption the eccentricity can 
depend on x as e = e + ejX 3 , with eo, ej being 
parameters, in order to encompass the usual oblate 
spheroid (ej = and < eo < 1), and a deformed 
oblate spheroid (eo > ej ^0). It is worth of calling 
attention that no dynamics is present if the second 
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exp(l/2Q a ) will be further approximated according 
to the schemes presented in the last Section: the 
first is to expand the exponential as a power series 
of ~Q/2 in the same way as shown in Eq. (jTTJ) ; the 
second is to expand the exponential using the col- 
location method. In order to provide a quantitative 
measure of the error between these approximations 
and the exact expressions for the corresponding ini- 
tial data, the L2 error defined by 




[K Q (x)-K a (x)]' 2 dx (15) 



will be evaluated in each case. 
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(b) 



FIG. 1: L2 error for the oblate spheroid (e^ = 0) in func- 
tion of the eccentricity e for J = 1, 2, 3, 4 (upper to lower 
curves). Similar plots are constructed in (b) considering 
the second scheme of approximation, for which the we 
have selected M — 6,8 and 10 collocation points start- 
ing from the upper curve (boxes). Note that the last two 
curves almost coincide indicating a rapid convergence of 
the method, whereas in (a) the convergence is evident 
only if J > 3. In all cases we have set N = 7. 



data family is chosen, since it represents an exact 
stationary solution of the RT equation. 

According to the previous discussion we choose 
the truncation order N from which the decomposi- 
tion Q a (uQ,x) — J2k=obk(u )Pk(x) is established. 
The values of each modal coefficient fefc(uo) will 
be determined using the standard Galerkin projec- 
tion procedure starting from Kq(x) ~ K a {x) = 



A exp(l/2 J2k=o bk(u)Pk(x)) that yields 



bk(uo) 



(2\n(Ao 1 K (x),P k (x)) 
(P k (x),P k (x)) 



(14) 



Henceforth, we set Aq = ko = 1 without loss of gen- 
erality. Once the modal coefficients are determined 




(a) 




FIG. 2: Z/2 error for the deformed oblate spheroid (e.2 = 
0.07) in function of the eccentricity eo for (a) J = 2, 3, 4 
(upper to lower curves). Similar plots are constructed 
in (b) considering the second scheme of approximation, 
for M — 6, 8 and 10 collocation points starting from the 
upper curve. Note that as shown in Fig. 1(b) the last two 
curves almost coincide indicating a rapid convergence of 
the method. In all cases we have set N = 7. 



Let us consider the oblate spheroid with ej = 
and < e < 1. We then fix the truncation 
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order to N = 7 and evaluate the error L 2 be- 
tween the exact and approximate expressions for 
the oblate spheroid taking into account both addi- 
tional schemes as previously described. In the first 
we choose J = 1,2,3 and 4 present in the expan- 
sion exp(Q /2) ~ J2k=o 2~ k Qa/&, whereas for the 
second scheme the expansion of exp((5 a /2) is con- 
structed with M = 6, 8 and 10 collocation points. 
The respective L 2 errors evaluated for successive val- 
ues of the eccentricity e are shown in Fig. 1. Two 
important aspect are worth to be mentioned. As ex- 
pected small eccentricities produce very small errors 
even for the quite modest choice J = 1, since these 
cases correspond to tiny departures from the spher- 
ically symmetric configuration described exactly by 
the Galerkin decomposition as &o(wo) = constant, 
bk(uo) = for k ^ 0. High eccentricities, otherwise, 
represent more severe deviations from the spheri- 
cal configuration, producing as a consequence an in- 
crease of the error. The second aspect concerns to 
the fast convergence achieved using both approaches 
as suggested by the plots of Figs. 1(a) and 1(b). Ac- 
cordingly, in the first case the convergence is clear if 
J > 3, and it is necessary that the number of collo- 
cation points is greater than, or N > 8 1 . Despite the 
acceptable errors attained with a modest truncation 
order (N = 7) in both cases, the use of the second 
scheme appears to be superior because the less alge- 
braic effort necessary to achieve the same results of 
the first case with J = 4. 

In Fig. 2 we have evaluated the L 2 error for the 
(deformed) oblate spheroid with variable eccentric- 
ity e — cq + e 2 x 2 , in which we have set e 2 — 0.07 
and again < e < 1. As it can be seem from the 
plots it is clear that the scheme using the collocation 
points is more efficient, in the sense of the need of 
smaller computational effort, and accurate than the 
first scheme. The initial data representing a black 
hole in motion is considered in Fig. 3, where J = 4 
and M = 8 collocation points which confirms once 
more the previous conclusion. 

In the second group of numerical experiments we 
shall verify the evolution of the error associated to 
the constraint equation ^j. The procedure consists 
into rewrite Eq. @ in the following way 
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FIG. 3: L2 error for the initial data identified as the 
boosted black hole (|13[) if the truncation order is TV = 7, 
and taking into account the collocation scheme (M = 8, 
lower curve) and the first scheme with J = 4 (upper 
curve) for the first scheme. 



Here X(u,x) is given by Eq. ©, but K(u,x) will 
be expressed according to the schemes of approxi- 
mation as outlined by Eqs. ifTTj) and (p~2|) . For in- 
stance, by assuming the expansion provided by the 
later expression we have 



Vfe=o / 

and, as a consequence, the constraint equation now 
reads 



CE(u, x) = X a A 2 



\qk{u)T k (x) 



/ 1 _ a ,2N.Efelo 1k{u)T' k {x) 



Efci *( M ) T fe( a 



- 1 = 0.(17) 



-1 = 0. 



(16) 



A similar expression is obtained if the approxima- 
tion expansion given by (jllj) is assumed. In any 
case the resulting form of the constraint equation is 
evaluated at each instant after the integration of the 
dynamical equation for the modal coefficients. As a 
matter of fact, we are interested in the rms error of 
the constraint or 



1 The convergence depends on the truncation order, but we 
have noticed that for a fixed truncation order N the most 
efficient decomposition is obtained if M = N+l collocation 
points is selected. If M > N + l there is not relevant 
improvement in the error as exemplified in Fig. 1(b). 




CE(u,x) 2 dx 



(18) 



Fig. 4 shows the evolution of this error if the initial 
data is the oblate spheroid (eo = 0.8, ej = 0) for 
both schemes of approximation. 



() 




(a) 




In Fig. 5 the convergence is illustrated by graphs 
corresponding to several truncation orders, namely, 
N = 5, 7,9, 11, and respectively M + 1 collocation 
points. As it can be seem the error decreases rapidly 
until reaching to the value considered zero up to our 
numerical precision, and the greater is the trunca- 
tion order less time is necessary to attain this value. 




FIG. 5: Evolution of the residual error the truncation 
orders N = 5,7, 9, 11 (upper to lower curves) and respec- 
tive N + 1 collocation points. The initial conditions cor- 
respond to the oblate spheroid with eccentricity eo = 0.8. 



(b) 



FIG. 4: (a) Evolution of the L2 error of the con- 
straint equation using the first and the second schemes 
of approximation (upper and lower curves, respectively). 
Here N = 7, J = 2 (see Eq. (JTTJl) and M = 8 collocation 
points, where it is clear the rapid convergence of the col- 
location expansion (|12[) . In (b) we present the influence 
of increasing the truncation orders N = 5, 7, 9, 11 with 
respective N + 1 collocation points. The initial condi- 
tions correspond to the oblate spheroid with eccentricity 
e = 0.8. 



The residual equation can provide another possi- 
ble way of verifying the convergence and accuracy 
of the numerical method as the truncation order is 
increased. In fact, as demonstrated by ??, the rms 
error of residual equation is identified as the upper 
bound of the rms error between the approximate and 
exact solutions. Then, in our last numerical experi- 
ment the evolution of the error 




Res(tt, x) 2 dx, 



where Res(u, x) is the residual equation (cf. Eq. ©) 
considering only the second scheme of approxima- 
tion (Eq. (|12|)) based on the collocation method. 



V. FINAL REMARKS 

The low dynamical system approach provided by 
the Galerkin method constitutes a valuable tool to 
deal with nonlinear problems as demonstrated not 
only here, but also in problem of turbulence of 
fluids 5]. Two aspects are worth of mentioning: 

• The transformation of, say, a partial differen- 
tial equation to a finite set of ordinary differ- 
ential equations. In general the greater is the 
set of these equations more accurate is the ap- 
proximate solution. 

• A reasonable accuracy can in fact be obtained 
with a relatively low dimensional dynamical 
system. 

As a matter of fact we have applied the Gakerkin 
method to a problem of gravitation, namely, the dy- 
namics of the Robinson- Trautman spacetimes that 
represents the simplest axisymmetric geometries en- 
dowed with gravitational waves. As we have de- 
scribed one of the crucial steps to implement the 
method is to project the residual equation into each 
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basis function, which is necessary to obtain the dy- 
namical system. Such projections are integrations 
in the spatial domain (in our case angular domain), 
but due to the presence of an exponential term, e~® , 
these integrations can not be performed as desirable. 
Then, we have engendered two schemes to overcome 
this difficulty (cf. Eqs. (HH) and ([12])). where the 
decomposition of the exponential with respect to a 
series of Chebyshev functions using the collocation 
method appeared to be quite satisfactory as demon- 
strated by several numerical experiments. Accord- 
ingly, even for a small truncation order, say N = 7 



(and N + 1 collocation points) to decompose e~® , 
the error associated to the constraint and residual 
equation are acceptable. The increase of truncation 
error to N = 11 - still a modest truncation order 
- provides a much better accuracy and fast conver- 
gence. 

Finally, the combination of Galerkin and colloca- 
tion methods can be very useful in studying other 
interesting nonlinear problems in Cosmology and 
Gravitation. 
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